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P^ i Abstract 

\—{ \ Newton's method for solving the matrix equation F{X) = AX — 

^^' XX AX = runs up against the fact that its zeros are not isolated. 

This is due to a symmetry of F by the action of the orthogonal group. 
C^^ I We show how differential-geometric techniques can be exploited to re- 

move this symmetry and obtain a "geometric" Newton algorithm that 
•^ ■ finds the zeros of F. The geometric Newton method does not suffer 

^" ! from the degeneracy issue that stands in the way of the original Newton 

f-? ' method. 

-(— > 
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^ ! 1 Introduction 

OO . Let A be a symmetric positive-definite n x n matrix and let p be a positive integer smaller 

than n. 

Oja's flow |Oja82[[Qja89| , in its averaged version 



o 



o, d 

g . ^X(t) = F{X{t)), (la) 

^ ' where F denotes the vector field 

j!^ ■ F: W'P -^ M"^P :X^AX- XX'^AX, (lb) 

C^ ■ 
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is well known for its principal component analysis properties. For all initial conditions Xq, 
the ordinary differential equation ([T]) has a unique solution curve t i— > X{t) on the interval 
[0,00) IYHM941 Th. 2.1], the limit X(cx3) = limi^oo^(*) exists, the convergence to X(oo) is 
exponential, and X(oo) is a zero of Oja's vector field F (jlbp |YHM94l Th. 3.1]. 
Observe that the zeros of F are the solutions X £ M"^^ of the matrix equation 

AX = XX^AX, 

which implies that the columns of AX are linear combinations of the columns of X. Letting 

col(y) = {Ya :aeW} 

denote the column space of y € R"^^, it follows that all zeros X of (jlbp satisfy Acol{X) C 
col(X), i.e., col(X) is an invariant subspace of A. Moreover, X(oo) is an orthonormal ma- 
trix (X (oo)X(oo) = /p), thus of full rank, whenever the initial condition X{0) is of full 
rank |YHM941 Prop. 3.1]. Assuming that X{0) has fuh rank, it holds that co1(X(cxd)) is the 
p-dimensional principal subspace of A if and only if col(X(0)) does not contain any direc- 
tion orthogonal to that subspace |YHM94| Th. 5.1]. The set of all initial conditions that 
do not satisfy this condition is a negligible set, i.e., Oja's flow asymptotically computes the 
p-dimensional principal subspace of A for almost all initial conditions. We also point out that 
Oja's flow induces a subspace flow, called the power flow |ASM08] . 

Because of these remarkable properties, Oja's flow has been and remains the subject 
of much attention, in its form ([T]) as well as in several modified versions; see, for exam- 
ple, IYHM941 lYan98l IDMV9"9l IMHM051 IJO06J . However, turning Oja's flow into a principal 
component algorithm implementable in a digital computer requires to discretize the flow in 
such a way that its good convergence properties are preserved. Since the solutions X{t) con- 
verge exponentially to their limit point X(oo), it follows that the sequence of equally spaced 
discrete-time samples {X{kT))k£N converges only Q-linearly |OR70] to X(oo). Therefore, nu- 
merical integration methods that try to compute accurately the solution trajectories of Oja's 
flow are not expected to converge faster than linearly. 

Nevertheless, if the ultimate goal is to compute the principal eigenspace of A, then it 
is tempting to try to accelerate the convergence when the iterates get close to the limit 



point X(oo), using techniques akin to those proposed in Hig99 IFK05J KQL+06 ILKLT06] . 



in order to obtain a superlinear algorithm. To this end, it is interesting to investigate how 
superlinearly convergent methods perform for finding a zero of Oja's vector field (jlbp . Since 
Newton's method can be thought of as the prototype superlinear algorithm to which all other 
superlinear algorithms relate, we propose to investigate the behavior of Newton's method 
applied to the problem of finding a zero of Oja's vector field (|lbp . 

A crucial hypothesis in the classical local convergence analysis of Newton's method (see, 
e.g., |DS83j ) is that the targeted zero is nondegenerate. As it turns out, the zeros of Oja's 
vector field F (jlbp are never isolated, because F displays a property of symmetry under 
the action of the orthogonal group Op on the set M"^^. Therefore, the classical superlinear 
convergence result of Newton's method is void in the case of F, and indeed our numerical 
experiments show that Newton's method on M"^^ for F performs poorly (see Figure [1]). 

In this paper, we propose a remedy to this difficulty, that consists in "quotienting out" the 
symmetry of F. Conceptually speaking, instead of performing a Newton iteration in M"^^, 
we perform a Newton iteration on a quotient space, namely M" ^/Op (defined in Section [3]). 
We exploit the Riemannian quotient manifold structure of M* ^/Op to formulate a Newton 



method on this set, following the framework developed in |AMS08j . The resulting Newton 
equation is a linear matrix equation that can be solved by various numerical approaches. 
It follows from the theory of Newton method on manifolds (see, e.g., [ADM^02| IAMS08J ). 
and from a careful analysis of the zeros of the vector field, that the obtained algorithm 
converges locally superlinearly to the set of orthonormal bases of invariant subspaces of A. 
Our numerical experiments show that the method behaves as expected. 

2 Plain Newton method for Oja's vector field 

We assume throughout the paper that j4 is a real symmetric positive-definite n x n matrix. 
For simplicity, we also assume that the eigenvalues of A satisfy 

Ai > • • • > A„, (2) 

i.e., all the eigenvalues of A are simple. Hence the p-dimensional invariant subspaces of A are 
isolated. 

Newton's method in M^^^ for Oja's vector field F ()lb|) consists of iterating the map 
X I— > X^ defined by solving 

BF{X)[Z] =AZ- ZX'^AX - XZ'^AX - XX'^AZ = -F{X) (3) 

X+ = X + Z. (4) 

The following proposition is a well-known characterization of the zeros of F. 



Proposition 2.1 (zeros of F ()lb|) ) Let X G IR"^^ be of full rank. Then the following two 
conditions are equivalent: 

1. F{X) =0, i.e., 

AX = XX'^AX. (5) 

2. col(X) is an invariant subspace of A and X is orthonormal (XX = I). 

Proof. 1^2. We have that AX = X{X'^AX), thus Aco\{X) C col(X). (More precisely, 
since X has full rank and A is positive-definite, it follows that X AX is invertible and 
thus Aco\{X) = col(X).) Moreover, multiplying © by X'^ on the left yields X'^ AX = 
X'^XX'^AX, which implies that X'^X = I since X'^AX is invertible. 

2 =^ 1. Since col(X) is an invariant subspace of A, there is a matrix M such that 
AX = XM. Multiply this equation on the left by X'^ to obtain that M = X'^AX and 
thus ©. D 

Hence, the set of all full-rank zeros of F is the finite union of the compact sets 

Si := {X G ]R"^P : col(X) = £i,X^X = 1} (6) 

where £i, . . . , £n are the p-dimensional invariant subspaces of A. (Note that N is finite; it is 
equal to (").) It is readily checked that Si = FjOp, where 

Op := {Q G WP : Q^Q = 1} 

is the orthogonal group of order p and Vi is an element of Si. It follows that the zeros of 
F are not isolated. Hence the Jacobian of F at the zeros of F is singular (i.e., the zeros are 



degenerate), and consequently, the classical result (see |DS83l Th. 5.2.1]) of local superlinear 
convergence of Newton's iteration to the nondegenerate zeros of F is void. This does not 
imply that Newton's method will fail, but there is a suspicion that it will behave poorly, and 
indeed, in Section (see Figured]), we report on numerical experiments showing that this is 
the case. 

3 Newton's method on ]R"^^/Op 

The degeneracy of the zeros of F is due to the following fundamental symmetry property: 

F{XQ) = F{X)Q, for all Q € Op. (7) 

In this section, we propose a geometric Newton method that performs well with functions 
F that satisfy d?]). This geometric Newton method evolves on the quotient space M"^^/Op, 
where this symmetry is removed. Then, in Section [H we will return to the specific case where 
F is Oja's vector field (jlbp and obtain a concrete numerical algorithm. 

The general idea for the geometric Newton method is first to define a vector field ^ on 
the manifold M" ^ /Op, whose zeros relate to those of F. The vector field ^ is specified in 
terms of its so-called horizontal lift in M" ^. This formulation requires us to introduce some 
concepts (vertical and horizontal spaces) borrowed from the theory of fiber bundles |KN63j , or 
more specifically from the theory of Riemannian submersions |Q'N83j . All the differential- 
geometric concepts used in this section are explained in |AMS08j . or in |Lee03| for what 
concerns Lie groups. 

The following notation will come useful. Let 

M^xp = {X £ W'P : det(X^X) ^ 0} (8) 

denote the set of all full-rank n x p matrices. Let 

symiB) = ^{B + B^) (9) 

and 

skew{B) = -{B - B^) (10) 

denote the terms of the decomposition of a square matrix B into a symmetric term and a 
skew-symmetric term. For X € M" ^, define 

pP ; ]R«xp ^ ^nxp . ^ ^ pp ^^) ^^j_ XiX^X)-^X^)Z (11) 

PI : Wixp ^ jgnxp . ^ ^ p|(^) ^ Xsym({X^Xr^X^Z) (12) 

P| : M"^P ^ M"^P : Z ^ P%{Z) = Xsk.e^{{X^ Xy^X'^ Z) (13) 

ph = pP + PI : Z h^ Z - Xskew((X^X)~^X^Z). (14) 

We have Z = P^{Z) + P^{Z) + P^{Z) for all Z G M"xp_ Observe that im(pP) = {Z e 
M.nxp . xTz = 0}, im(P|.) = XSsymip), im(Pl) = X5skew(p), where 

Ssymip) = {Se M^^P ■.S^ = S} 



p 



denotes the set of all symmetric matrices of order p and 

5skcw(p) = {^ e M^^*' : n^ = -n} 

is the set of all skew- symmetric (or antisymmetric) matrices of order p. (The letter "p" stands 
for "perpendicular" , "s" for symmetric, "a" for antisymmetric, and the notation P^ will make 
sense in a moment.) 

In M" ^, we define an equivalence relation "~" where X ~ y if and only if there exists a 
Q G Op such that Y = XQ. The equivalence class of X E M" ^ is thus 

[X] =XOp = {XQ : Q G Op}. (15) 

We let 

denote the quotient of M" ^ by this equivalence relation, i.e., the elements of the set M" ^ /O 
are the equivalence classes of the form (fT5]) . X G M" ^ . We let 

TT : M^P ^ M^VOp (16) 

denote the quotient map that sends X G M* ^ to its equivalence class [X] viewed as an element 
of M" ^ /Op. The set M* ^ is termed the total space of the quotient M" ^ /Op- Note that a 
point 7r(X) G M" '^/Op can be numerically represented by any element of its equivalence class 
[X]=^-i(7r(X)). 

Since M" ^ is an open subset of R"^^, it is naturally an open submanifold of the linear 
manifold M"^^. Moreover, it can be shown that the map 

^ : 0„ X R^^ ^ Mr^ : (Q, ^) ^ ^Q 

is a free and proper Lie group action on the manifold M" ^ . Therefore, by the quotient man- 
ifold theorem (see, e.g., |Lee03t Th. 9.16]), the orbit space M* ^ /Op is a quotient manifold. 
In other words, the set M" ^ /Op is turned into a manifold by endowing it with the unique 
differentiable structure that makes the quotient map vr a submersion. It comes as a conse- 
quence that each equivalence class [X] = tt~^{'k{X)), X G M^^^, is an embedded submanifold 
of M" ^ . We term vertical space at X G M" ^ the tangent space Vx to [X] at X, i.e., 

Vx = Tx[X] = X TjOp = X5skew(p). 

Observe that im(P^) = Vx- 

Next we define horizontal spaces fix, which must satisfy the condition that M"^^ is the 
internal direct sum of the vertical and horizontal spaces. We choose 

Hx = im(Pl) = {XS + X^K : S^ = S} (17) 

where P^ is as in (J14p and X± G ]R"'^("~p) denotes any orthonormal matrix such that X-^X_|_ = 
0. (It would be sufficient to state that X± has full rank, but it does not hurt to assume that 
it is orthonormal. We do not use X^ in the numerical algorithms.) 

As an aside, we point out that the horizontal space ()17p is the orthogonal complement of 
the vertical space with respect to the noncanonical metric g defined by 

g^iZi, Z2) = tr(Zf PP (Z2) + Zf X(X^X)-2x^Z2). 



Indeed, we have {Z G M">^p : g^(VF, Z) = for all VF E Vx} = {XS + X^K iS^ = S}. The 
reason for not choosing the canonical Riemannian metric on M" ^ is that it yields a more 
intricate formula for the horizontal space and for the projection onto it. 

Consider a function F : R"^p -^ W^^p that satisfies the symmetry property ([7]). Define a 
horizontal vector field ^ by 

^j, := P^{F{X)). 

It can be checked that this horizontal vector field ^ satisfies D7r(X)[^jf] = T)tt(XQ)[(^xq] 
for all Q G Op, and thus ^ is the horizontal lift of a vector field ^ on the quotient manifold 
M* ^ /Op. In the remainder of this section, we formulate a geometric Newton method for 
finding a zero of the vector field ^. 

For finding a zero of a vector field ^ on an abstract manifold M endowed with an affine 
connection V and a retraction R, we consider the geometric Newton method in the form 
proposed by Shub [ADM"'"02l §5] (or see |AMS08( §6.1]). The method consists of iterating 
the mapping that sends x G Al to x+ G Al obtained by solving 

hix)[nx] = -Cx, Vx G T^M, (18a) 

x+ = R^{r,^). (18b) 

Here J^(x) denotes the Riemannian Jacobian J({x) : TxM — *• TxM. : Cx "-^ Jc(2^)[Cx] = V^^^. 
For the case where A4 is the quotient M" ^/Op, we choose the affine connection V defined 

by 

(Vr,.,x)e)^:=^xDe(^)[^x], (19) 

where the over line denotes the horizontal lift. It can be shown that the right-hand side is 
indeed a horizontal lift and that the definition of V satisfies all the properties of an affine 
connection. With this choice for V, and with a simple choice for the retraction R (see |AMS08( 
§4.1.2] for the relevant theory), the geometric Newton method on the quotient manifold 
ffi* ^/Op becomes the iteration that maps vr(X) to 7r(X_|_) by solving 

PkB'i{X)[rix] = -Ix, rix^nx (20a) 

X+ = X + r)x. (20b) 



with P as in ()14p and 7ix as in (I17p . Note that, in spite of possibly unfamiliar notation, ()20p 
only involves basic calculus of functions between matrix spaces. 

The local superlinear convergence result for the geometric Newton method (I20p follows 
directly from the local convergence result of the general geometric Newton method (jlSp . 
see |AMS08t §6.3]. It states that the iteration converges locally quadratically to the nonde- 
generate zeros of ^. Since we have "quotiented out" the symmetry of F by the action of Op 
to obtain ^, it is now reasonable to hope that the zeros of ^ are nondegenerate. 

Note that the proposed geometric Newton method differs from the "canonical" Rieman- 
nian Newton algorithm |Smi94j on the quotient Mil ^/Op, because the affine connection chosen 
is not the Riemannian connection on M^ ^/Op endowed with the Riemannian metric inherited 
from the canonical metric in M" ^ . However, the property of local superlinear convergence 
to the nondegenerate zeros still holds (see |ADM"'"02J or |AMS08j ). 



4 A geometric Newton method for Oja's vector field 

In this section, we apply the geometric Newton method on M" ^ /Op, given in (|2U|) . to the 
case where the tangent vector field ^ on M" ^/Op is defined by the horizontal lift 

l^ := P^^{F{X)\ (21) 

with P^ as in (J14p and F as in (jlbp . The resulting Newton iteration is formulated in 
Algorithm m Recah the definitions of M"''^ ([8]), P'^ dH]), skew (ITO]). Hx (IT1 



Algorithm 1 Geometric Newton for Oja's vector field 



Require: Symmetric positive-definite n x n matrix A; positive integer p < n. 
Input: Initial iterate Xq E M" ^, i.e., Xq is a real n x p matrix with full rank. 
Output: Sequence of iterates {Xk) C M" ^. 

1: for /c = 0,1,2, ... do 

2: Solve the linear system of equations (we drop the subscript k for convenience) 

PxiAZ - ZX'^AX - XZ^AX - XX'^AZ - Zskew{{X'^X)-'^X'^AX) 
- Xskew{-{X^Xy^{X'^Z + Z^X){X^ Xy^X^AX + {X^X)~^{Z'^AX + X^AZ))) 

= - {AX - XX^AX - Xs'kew{{X'^X)-^X'^AX)) (22) 

for the unknown Z G Tix- 
3: Set 



^k+l — ^k + Z. 



4: end for 



Observe that Algorithm [T] is stated as an iteration in the total space M" ^ of the quotient 
M* ^ /Op. Formally, the sequence of iterates of the Newton method on M" ^ /Op, for an 
initial point 7r(Xo) € M* ^ /Op, is given by (vr(Xfc))fcgp^, where vr is the quotient map (fTOj) and 
(Xfc)fc£N is the sequence of iterates generated by Algorithm [TJ 

We point out that (j22p is merely a linear system of equations. It can be solved by means 
of iterative solvers that can handle linear systems in operator form. Moreover, these solvers 
can be stopped early to avoid unnecessary computational effort when the iterate X^ is still 
far away from the solution. Guidelines for stopping the linear system solver can be found, 
e.g., in |EW96j . In our numerical experiments (Section [5]) we have used Matlab's GMRES 
solver. 

Algorithm [1] converges locally quadratically to the nondegenerate zeros of ^. We first 
characterize the zeros of ^, then we show that they are all nondegenerate under the assump- 
tion ^. 

First note that ^^(x) = if and only if P^^{F{X)) = 0, where X € M"""^. Under the 



assumption that X e M* , the fohowing statements are equivalent: 

Pl(F(X))=0, 

F(X)Gim(Pl), 

F{X) = Xn for some = -0^, 

AX - XX'^AX = Xn for some Q = -fl^, 

Acol{X) C col(X) and {X'^Xy'^X'^AX - X'^AX is skew-symmetric, 

^col(X) C col(X) and {X'^Xy'^X'^AX + X'^ AX {X'^ X)"^ = 2X'^AX. 

We thus have an equation of the form YB + BY = 2B, where B := X^ AX is symmetric 
positive definite, hence its eigenvalues f3i, i = 1, . . . ,p, are all real and positive. The Sylvester 
operator Y i— > YB + BY is a linear operator whose eigenvalues are (3i + f3j, i = 1, . . . ,p, 
j = 1, . . . ,p |Gan591 Ch. VI]. All these eigenvalues are real and positive, thus nonzero, hence 
the operator is nonsingular, from which it follows that the equation YB + BY = 2B has one 
and only one solution Y. It is readily checked that this unique solution is y = /. We have 
thus shown that XX = I. The result is summarized in the following proposition. 



pnxp 
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Proposition 4.1 Let P^ be as in (fTil) and let F be Oja's vector field (fTb|) . Then X G 

a zero of the projected Oja vector field P^{F) if and only if Acol{X) C col(X) and X'^X = /. 

In other words, the full-rank zeros of the projected Oja vector field P (F) are the full-rank 
zeros of Oja's vector field F. This means that we do not lose information by choosing to 
search for the zeros of ^ (I2ip instead of F (llbh . 

It remains to show that the zeros of ^ are nondegenerate. Let 7r(X^,) be a zero of ^, 
which means that AX^ = X^XjAX^ and Xj^X^, = /. The task is to show that the Jacobian 
operator J^(7r(X*)), or equivalently its lifted counterpart 

J^(X,) : im(P^) - im(P^) : Z ^ P^ (d(P^F))(X.)[Z]) , (23) 

is nonsingular. To this end, consider the operator 

J : M«xp ^ j^nxp . ^ ^ D(P^(F))(X,)[Z]. (24) 

Note that J^{X^:) is the restriction of J to im(P'^). Consider the decomposition W^^p = 
im(PP) im(P'') © im(P'') and recall that im(P'') = im(PP) © im(P''). We show that the 
corresponding block decomposition of J is as follows: 



* 





0" 


? 


* 
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? 






where "*" denotes nonsingular operators. It then directly follows that the upper two-by-two 
block of J, which corresponds to J^(X*), is nonsingular. 

We show that the 11 block (i.e., the "pp" block) is nonsingular. This block is the operator 
from im(P^ ) to im(P^ ) given by 

Z ^ P^B{P''{F)){X,)[Z]=PP^AZ-P^^ZXJAX,. 



(In obtaining this result, we have used the relations Xj^X = I, skew[Xj AX^:) = 0, P^ X^ = 
0.) In view of the hypothesis that the eigenvalues of A are all simple, this operator is known 
to be nonsingular; see, e.g., |LE02tEEVM04j . 

We show that the 22 block (i.e., the "ss" block) is nonsingular. This block is the operator 
from im(P|- ) to im(P|' ) 

X,S ^ P!^(d{P^{F)){X,)[X,S])=-X,{SXJAX,+xJAX,S). 



(We have used the relation AX^, = X^X'^AX^ to obtain this expression.) In view of the 
previous discussion on the Sylvester operator, this operator is nonsingular. 

This completes the proof that the zeros of ^ are nondegenerate. Consequently, for all Xq 
sufficiently close to some Si ^ , the sequence (X^ ) generated by Algorithm [T] is such that 
XkOp converges quadratically to Si. Recall that Si is the set of all orthonormal matrices 
whose column space is the ith invariant subspace of A. 

5 Numerical experiments 

In this section, we report on numerical experiments for both the plain Newton and the ge- 
ometric Newton method, derived in Section 2 and Section 4, respectively. The experiments 
were run using Matlab. The machine epsilon is approximately 2 • 10~^^. 

As mentioned in Section 1, the plain Newton method performs poorly due to the fact that 
the zeros of the cost function F are not isolated. To illustrate this, we consider a symmetric 
positive definite matrix A G M^^^ with uniformly distributed eigenvalues in the interval [0, 1] 
and 100 different initial iterates Xq G M^^^. Each of the initial iterates is computed as 

Xo = X^ + IQ-^E , 

where X^, is such that F{X^) = and £^ is a 6 x 3 matrix with random entries, chosen from 
a normal distribution with zero mean and unit variance. The simulation is stopped after 50 
iterations. One representative example is given in Figure [H In Figure [1] (a), the norm of 
F{X) is given for each iteration step. Close to the solution, the system matrix gets singular 
and the algorithm deviates from the optimal point. In Figure [T](b), we present the evolution 
of the norms of the three components K, XQ, and XS of the update vector Z, 

z = x±K + xn + xs , 

where X± is the orthogonal complement of X, il is a skew-symmetric matrix and S" is a 
symmetric matrix. We see that, even when the K and XS component are very small, Xil. is 
quite large. This concords with the fact that the kernel of the Hessian at a stationary point 

X is {xn : n'^ = -n}. 

Next, we study the geometric Newton method, derived in Section 4. Again, we consider 
n = 6, p = 3 and a symmetric positive definite matrix A G M^^^ with uniformly distributed 
eigenvalues in the interval [0, 1]. We perform lO'^ experiments with a single matrix A but with 
different initial iterates Xq G K^^^ with random entries, chosen from a normal distribution 
with zero mean and unit variance. In Figure [21 we show the number of runs that converged 
to each of the eigenspaces of A. The dominant eigenspace is marked by "123". In general. 





(a)F 



(b) K, XS, Xn 



Figure 1: Plain Newton method 

"ijA;" stands for the eigenspace spanned by the i , j , and A;^ eigenvectors!^. Let W be the 
matrix of all eigenvectors. We declare that the algorithm has converged to eigenspace "ij/c" 
when the norms of the i^\ j^\ and k columns of the matrix X^W are all greater than 10"^'^ 
after 50 iterations and the norms of the rest of the columns are smaller than 10"^'^. It appears 
that the basin of attraction of the dominant eigenspace is the largest. In our experiment, all 
the runs have converged to one of the 20 possible eigenspaces. In general, there may be cases 
where the algorithm does not converge to any of the eigenspaces. This may occur when the 
initial iterate Xq is very close to the boundary of one of the basins of attraction. However, 
these cases are rare. Finally, the superlinear convergence rate of the algorithm is illustrated 
in Figure [3l 

6 Conclusion 

We have investigated the use of Newton's method to compute superlinearly the zeros of Oja's 
vector field (|lbp . Due to a symmetry in the vector field by the action of the orthogonal group, 
its zeros are never isolated, which causes the plain Newton method to behave poorly. We 
have proposed a remedy that consists in "quotienting out" the symmetry. This led to the 
formulation of a geometric Newton algorithm that seeks the zeros of a projection of Oja's 
vector field. We have shown that the zeros of the projected vector field are the same as the 
zeros of the original vector field. Moreover, these zeros are nondegenerate. This means that by 
quotienting out the action of the orthogonal group, we have removed just enough symmetry 
to make the zeros nondegenerate. In view of the nondegeneracy property, it follows directly 
from the convergence theory of the abstract geometric Newton method that the resulting 
algorithm converges locally superlinearly to the zeros of Oja's vector field. 

Invariant subspace computation has been and still is a very active area of research. 
As a method for invariant subspace computation, it is doubtful that the proposed algo- 



^For convenience, we consider an eigenvalue decomposition, where the eigenvalues of A axe put in descending 
order. 
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Figure 2: Geometric Newton method 

rithm can outperform the state-of-the-art methods. In particular, the Grassmann-based ap- 
proach |EAS981 ILE021 IAMS04J , that can be thought of as quotienting out the action of the 
whole general linear group instead of the orthogonal group, leads to a Newton equation that 
lives in a smaller subspace of MJ^^p and that can be solved in fewer flops. When n ^ p, how- 
ever, the number of operations to compute the iterates is of the same order. Moreover, the 
Grassmann-based algorithm does not possess the remarkable feature of naturally converging 
towards orthonormal matrices, i.e., without having to resort to orthogonalization steps such 
as Gram-Schmidt. 

The problem of computing the zeros of Oja's vector field was also an occasion for intro- 
ducing the quotient manifold M" ^/Op, that seems to have received little attention in the 
literature, in contrast to the more famous Grassmann and Stiefel manifolds. In later work, 
we will further analyze the geometry of this quotient manifold, which was just touched upon 
in this paper, and we will show how it can be used in the context of low-rank approximation 
problems. 
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